Chapter 19 Sympy
From book Python Programming (Problem solving, Packages and Libraries) published by McGraw Hill Education (India) Private limited.
By:

  • Anurag Gupta
  • G. P. Biswas

Note the following:-

  1. This html document is meant as an accompaniment to Chapter 19 Sympy .
  2. The document contains scripts executed on IDLE as well as on Jupyter notebook.
  3. The scripts executed on Jupyter can be directly copied and run into a Jupyter notebook or some other IDE (Like Pycharm or Eclipse with PyDev or Visual studio).
  4. However the scripts on IDLE also contain the >>> symbol and therefore cannot be directly executed. If you want to execute them on IDLE or Jupyter, you need to manually remove the >>> symbol.
  5. Wherever needed some background material from the book is also included to help you better understand the scripts
  6. The topic numbers given on each paragraph match the topic numbers of the book, so you can easily identify the topics and corresponding scripts.
  7. In some of the scripts, the file paths give are that of the author's computer. You need to replace them with file paths of your own computer.
  8. At some places, to improve readability, page numbers of the book are indicated in green font like:- See Page 181 of the book
  9. This document was first created as a Jupyter Notebook as combination of Markdown and code cells (extension .ipynb) and then downloaded as html. If someone wants to "modify" or "extend' this document, you may ask for the original .ipynb file by sending me an e-mail at:- 999.anuraggupta@gmail.com

19.2. Basics of SymPy 1
19.2.1. The “symbols()” function
In SymPy there is a function “symbols()” by which you can declare a variable to be a symbol. A variable “declared to be a symbol” using the symbols() function is like any other Python variable and, therefore, must be “assigned” before they can be used. This will be clear from the following code:
This script is available on page 481 of the book

In [1]:
from sympy import symbols
x,y,z = symbols('x y z')
a_expr = (x + y)*(y + z)
print(a_expr)

(x + y)*(y + z)

To help the programmer write the code faster, the symbols() function supports what is called the “range index”. The “range index” is indicated by a colon, i.e., (:). Further the “type of range” is determined by the “type of character” to the right of the colon. This will be clear from following code:

In [2]:
print(symbols('x:5'))
print(symbols('x10:15'))

(x0, x1, x2, x3, x4)
(x10, x11, x12, x13, x14)

The SymPy library also has expand and factor functions to expand/ factorize the expression as shown in the following code:

In [3]:
from sympy import expand, factor
exp1 = (x + y)*(x + y)*(y + z)
exp2 = expand(exp1)
print('exp1 on expansion->', exp2)
exp3 = factor(exp2)
print('exp2 on factorization->', exp3)

exp1 on expansion-> x**2*y + x**2*z + 2*x*y**2 + 2*x*y*z + y**3 + y**2*z
exp2 on factorization-> (x + y)**2*(y + z)

19.2.2. Importing symbols from module sympy.abc
The SymPy library has a module abc and one can directly import symbols from this module also.
This becomes clear if you look at the docstring of the abc module, which is also shown below (On Jupyter):

In [4]:
import sympy
?sympy.abc

Object `sympy.abc` not found.

19.2.4. Equality testing in SymPy using “$==$”
Just like Python, the symbol “$==$” is used for equality testing in SymPy. But many a times $a==b$ in SymPy may give unpredictable results.
Therefore, it is better to test whether $a-b ==0$.
This is done by a function called simplify. This will be clear from the following example:
This script is available on page 482 of the book

In [5]:
from sympy import *
a = (x + y)*(x - y)
b = x**2 - y**2
print(' Is a == b?', a == b)
print('is simplify(a - b) ==0?', simplify(a - b) ==0)

 Is a == b? False
is simplify(a - b) ==0? True

19.2.6. Using operators on combination of SymPy objects and Python objects
When two objects are linked by an operator and they are not of the same type, then some implicit type casting takes place. So, you may have two SymPy objects, two Python objects or one of each.
The problem arises when you want to write a formula x + ¾. Here, ¾ will evaluate to 0.75. SymPy also has an Integer class which converts a Python integer to a SymPy integer.
The way to write this is as follows:
This script is available on page 483 of the book

In [6]:
from sympy import *
a = x + 3/4
print(a) # Gives x + 0.75
b = x + Rational(3,4)
print(b)# Gives x + 3/4
c = x + Integer(3)/ Integer(4)
print(c) # Gives x + ¾

x + 0.75
x + 3/4
x + 3/4

19.3. Basics of SymPy 2
19.3.1. Substitution in a SymPy expression
In SymPy it is possible to replace a symbol with another symbol or SymPy object or Python object.
The following example will clarify the concept:
This script is available on page 484 of the book

In [7]:
from sympy import *
x, y, z = symbols('x, y, z')
a, b, c = symbols('a b c')
m = x**2 + y*4 + sin(z)
n = m.subs([(x, a), (y, 3), (z, 2*c)])
print(n) # Gives a**2 + sin(2*c) + 12

a**2 + sin(2*c) + 12

19.3.2. Convert Python strings to SymPy expression and evaluating it (Functions sympify() and evalf())
SymPy has a function sympify() (not simplify). sympify() can be used to convert strings to expressions and a function evalf() to get its value. Following example clarifies this:
This script is available on page 484 of the book

In [8]:
from sympy import *
#x, y, z = symbols('x y z')
a = 'x**2 + log(y) + cos(z)'
b = sympify(a)
print(b) # a -> x**2 + log(y) + cos(z). Note in sympy log is natural log ie ln
c = b.subs([(x,2), (y, 100), (z, 0)])
print(c) # Prints-> log(100) + 5
d = c.evalf()
print(d) # Output-> 9.60517018598809

x**2 + log(y) + cos(z)
log(100) + 5
9.60517018598809

sympify() can also be used to convert numbers of Python data type to SymPy data type. For instance, the number 1 in Python is an integer, while in SymPy it may be a rational. ½ in Python is 0.5 while Rational(1,2) or sympify(1)/2 in SymPy is ½. This is clear from the following example:
This script is available on page 484 of the book

In [9]:
from sympy import *
print('1/2 ->', 1/2)
print('Rational(1,2)->', Rational(1,2))
print('sympify(1)/2->', sympify(1)/2)

1/2 -> 0.5
Rational(1,2)-> 1/2
sympify(1)/2-> 1/2

19.3.3. Singleton class in SymPy
In SymPy there is a class called SingletonRegistry. Here we would not go into details but touch the essential factors:

  • Common mathematical constants are represented by Singleton classes.
  • This singleton class can be represented by capital letter S. So, S is an “instance ” of SingletonRegistry class.
  • These constants are directly available in the SymPy namespace. So, if you do an import sympy*, then S becomes available and through S, the constants also become available. Therefore, constants, such as “pi”,”e” can be accessed as S.pi, S.e, etc.
  • Infact, S can also act as a shortcut for sympify function. Remember that in SymPy if you want to get ½ as ½ and not 0.5, then use Rational(1,2). You can also use S(1)/2.

The following script shows that S(1)/2 is same as sympify(1)/2
This script is available on page 485 of the book

In [10]:
import sympy
from sympy import *
print('using S(1)/2->', S(1)/2)
print('Using sympify(1)/2->', sympify(1)/2)

using S(1)/2-> 1/2
Using sympify(1)/2-> 1/2

Table 19.2: Some common singletons of SymPy
Table 19.2 in the book shows some common siongletons of Sympy. The following example code from the book demonstrates their usage:-
This script is available on page 485 of the book

In [11]:
from sympy import *
print(1/2) # gives 0.5
print(Rational(1,2)) # gives 1/2
print(S(1)/2) # gives 1/2
print(pi, S.Pi) # Both pi and S.Pi give pi
print(E, S.Exp1) # Both E and S.Exp1 give E
print(oo, S.Infinity)# Both oo and S.Infinity give oo

0.5
1/2
1/2
pi pi
E E
oo oo

You can also see the complete list of all singletons (and other attributes) available in S by using the dir(S) function as shown below:
This script is available on page 486 of the book

In [12]:
from sympy import *
print(dir(S))

['Catalan', 'ComplexInfinity', 'Complexes', 'EulerGamma', 'Exp1', 'GoldenRatio', 'Half', 'IdentityFunction', 'ImaginaryUnit', 'Infinity', 'NaN', 'Naturals0', 'NegativeInfinity', 'NegativeOne', 'One', 'Pi', 'Reals', 'Zero', '__call__', '__class__', '__delattr__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattr__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__sizeof__', '__slots__', '__str__', '__subclasshook__', '_classes_to_install', 'false', 'register', 'true']

In Python, all objects have a method called mro(). Using this method you can know the parentage of a class in multiple inheritances. It returns a list of all the parent classes of the object. The order is by depth meaning that the shallowest class is listed first followed successively by each parent. For instance, the following code shows the mro() of S(1):
This script is available on page 486 of the book

In [13]:
from sympy import *
print('S(1).mro()->', type(S(1)).mro())

S(1).mro()-> [<class 'sympy.core.numbers.One'>, <class 'sympy.core.numbers.IntegerConstant'>, <class 'sympy.core.numbers.Integer'>, <class 'sympy.core.numbers.Rational'>, <class 'sympy.core.numbers.Number'>, <class 'sympy.core.expr.AtomicExpr'>, <class 'sympy.core.basic.Atom'>, <class 'sympy.core.expr.Expr'>, <class 'sympy.core.basic.Basic'>, <class 'sympy.core.evalf.EvalfMixin'>, <class 'object'>]

19.3.4. Functions in SymPy As per the SymPy documentation , there are three types of functions available in SymPy:

  • Defined functions (in the sense that they can be evaluated) such as exp or sin; they have a name and a body: f = exp
  • Undefined functions which have a name but no body. Undefined functions can be defined using a Function class as follows: f = Function('f'). (the result will be a Function instance).
  • Anonymous function (or lambda function) which have a body (defined with dummy variables) but have no name. Example: For 1 variable→ f = Lambda(x, exp(x)*x), for 2 variables→ f = Lambda((x, y), exp(x)*y).

The fourth type of functions are composites, such as (sin + cos)(x); these work in SymPy core, but are not yet part of SymPy.”
You can see how the undefined function is implemented in SymPy in the following code:
This script is available on page 487 of the book

In [14]:
from sympy import *
x,y = symbols('x y')
f = Function('f')
g = Function('g')(y)
h = Function('h')(x,y)
print('f->', f, 'g->', g, 'h->', h)# f-> f g-> g(y) h-> h(x, y)
j = f(x)
print(f(x), f(y))# Prints f(x) f(y)

f-> f g-> g(y) h-> h(x, y)
f(x) f(y)
**Advanced discussion on undefined functions: ** In the following code example, `f` is an undefined function but `g` is not. This can be confirmed from the type of `f` and `g`. Rather `g` is of type `AppliedUndef`. Come to think of it, this is somewhat like a function say sin or cos or log. By itself sin or cos are UndefinedFunction, but `sin(x)` or `cos(x)` are what would be `AppliedUndef` functions.

This script is available on page 487 of the book

In [15]:
from sympy import *
f = Function('f')
# type f is UndefinedFunction
print('type f->', type(f))
x = symbols('x')
# but type g is g
g = Function('g')(x)
print('type g->', type(g))
# lets see the class hierarchy of f and g
print('mro of f->', type(f).__mro__)
print('mro of g->', type(g).__mro__)

type f-> <class 'sympy.core.function.UndefinedFunction'>
type g-> g
mro of f-> (<class 'sympy.core.function.UndefinedFunction'>, <class 'sympy.core.function.FunctionClass'>, <class 'sympy.core.assumptions.ManagedProperties'>, <class 'sympy.core.core.BasicMeta'>, <class 'type'>, <class 'object'>)
mro of g-> (g, AppliedUndef, Function, Application, <class 'sympy.core.expr.Expr'>, <class 'sympy.core.basic.Basic'>, <class 'sympy.core.evalf.EvalfMixin'>, <class 'object'>)

19.3.5. Lambda class in SymPy
SymPy Lambda class is different from lambda function in Python. A good explanation of Lambda class and its use is given in SymPy docs :
“Lambda(x, expr) represents a lambda function similar to Python's lambda x: expr. A function of several variables is written as Lambda((x, y, ...), expr).” The docstring of the Lambda class gives a number of examples on how to use this Lambda class. On Jupyter, the docstring can be accessed using (?). A modified version of the docstring of SymPy’s Lambda class is shown below and some relevant lines of code are also explained:
This script is available on page 488 of the book

In [16]:
from sympy import Lambda
?Lambda

# OUTPUT (Truncated and modified)
Init signature: Lambda(variables, expr)
Docstring:     
Lambda(x, expr) represents a lambda function similar to Python's 'lambda x: expr'. A function of several variables is written as Lambda((x, y, ...), expr).

A simple example:
>>> from sympy import Lambda
>>> from sympy.abc import x
>>> f = Lambda(x, x**2)
>>> f(4)
16

For multivariate functions, use:
>>> from sympy.abc import y, z, t
>>> f2 = Lambda((x, y, z, t), x + y**z + t**z)
>>> f2(1, 2, 3, 4)
73

A handy shortcut for lots of arguments:

>>> p = x, y, z
>>> f = Lambda(p, x + y*z)
>>> f(*p)
x + y*z

19.4. Sets in SymPy
There is a detailed discussion in the book on the topic of Sets starting from page 488
19.4.1. FiniteSet
The FiniteSet can “contain” numbers. The following code shows how FiniteSet() works:
This script is available on page 490 of the book

In [17]:
from sympy import *
f1 = FiniteSet(1,2,3,4)  # Create FiniteSet by giving it some items
print(f1)
a_list = [2,3,4,5]
f2 = FiniteSet(*a_list)  # Create FiniteSet from a list
print(f2)
a_tup = (5,6,7,8)
f3 = FiniteSet(*a_tup)  # Create FiniteSet from a tuple
print(f3)
is_member = 5in f3 #Check whether an object member or not
print(is_member) # True

{1, 2, 3, 4}
{2, 3, 4, 5}
{5, 6, 7, 8}
True

The following code shows how to iterate over a FiniteSet:
This script is available on pages 490-491 of the book

In [18]:
from sympy import *
s1 = FiniteSet('A', 'B', 'C')
check_s1 = s1.is_iterable
print(check_s1)
if check_s1: #Should check whether iterable
    memb_s1 = iter(s1)
    print(memb_s1)
    len_s1 = len(s1)# Get number of members
    for r in range(len_s1):
        a_memb = next(memb_s1)# all iterables implement next()
        print(a_memb)

True
<tuple_iterator object at 0x05FFC7F0>
A
B
C

19.4.2. Interval
The signature of Interval class, which is an extension of Python set data type is:

Interval(start, end, left_open=False, right_open=False)

So by default both the left and right intervals are not open, that is, closed. Only real end points are supported (not complex). It must also be kept in mind that Interval(a, b) with a > b will return the empty set. The following code clarifies the concept:
This script is available on page 491 of the book

In [19]:
from sympy import *
print(Interval(0,1))# Equivalent to [0,1]
print(Interval(0,1, False, False)) # Equivalent to [0,1]
print(Interval(0,1, False, True))# Equivalent to [0,1). Interval.Ropen(0, 1)
print(Interval(0,1, True, True)) # Equivalent to (0,1). Interval.open(0, 1)
print(Interval(0,1, True, False))# Equivalent to (0,1]. Interval.Lopen(0, 1)
print(Interval(1,0))# EmptySet() 

Interval(0, 1)
Interval(0, 1)
Interval.Ropen(0, 1)
Interval.open(0, 1)
Interval.Lopen(0, 1)
EmptySet()

19.4.3. EmptySet
Following code shows how EmptySet() works. Intersection with an EmptySet() is always an EmptySet().
This script is available on page 491 of the book

In [20]:
from sympy import *
print(S.EmptySet) # Output is EmptySet()
f1 = FiniteSet(1,2,3)
print(f1.intersect(S.EmptySet))# Intersection with EmptySet() is EmptySet()

EmptySet()
EmptySet()

19.4.4. Intersection
Intersection(Set_A, Set_B) will return the intersection set of the two sets. You can also use the method intersect as set_A.intersect(set_B).
This is clear from the following code:
This script is available on page 491 of the book

In [21]:
from sympy import *
s1 = Intersection(Interval(2,4), Interval(3,5))# Use Intersection() class
print(s1)
s2 = Interval(2,4).intersect(Interval(3,5))# use intersect() method
print(s2)

Interval(3, 4)
Interval(3, 4)

19.4.5. Union

Short note on various printing options in SymPy: SymPy provides a number of “printers”. For details seethe link . To get the best pretty printing, you should use the init_printing() function. This function on its own enables the best printer available in your environment. The following code will initiate pretty printing on your machine:But this code works only if the terminal supports Unicode. If it doesn’t, you may use the ASCII pretty printer through pprint(). 

from sympy import *
init_printing()

Union(set_A, set_B) returns the union of “set_A” and “set_B”. As a shortcut one can use the “+” operator for “union” of two sets.
This script is available on page 492 of the book

In [22]:
from sympy import *
s1 = Union(Interval(0,1), Interval(2,3))
pprint(s1) # output is [0, 1] ∪ [2, 3]
s2 = Union(Interval(0,1), Interval(1,2))
pprint(s2) #Output is [0, 2]
s3 = Union(Interval(0,1, True, True), Interval(1,2, True, True))
pprint(s3)#Output is (0, 1) ∪ (1, 2)
s4 = Interval(0,1, True, True).union(FiniteSet(1,2))
pprint(s4)# Output is (0, 1] ∪ {2}

[0, 1] ∪ [2, 3]
[0, 2]
(0, 1) ∪ (1, 2)
(0, 1] ∪ {2}

19.4.6. ConditionSet
ConditionSet is a set that satisfies a given condition. The signature of ConditionSet on Jupyter is shown in the following example:
This script is available on page 492 of the book

from sympy import *
?ConditionSet

The output on Jupyter is somewhat like this:-

Init signature: ConditionSet(sym, condition, base_set)
Docstring:     
Set of elements which satisfies a given condition.
{x | condition(x) is True for x in S}

The use of ConditionSet will become clear from the code below:-
This script is available on page 493 of the book

In [23]:
from sympy import *
from sympy.abc import x
s1 = ConditionSet(sym = x,   # symbol is x
                  condition = Eq(x**2, 4 ),# equation is x **2 = 4
                  base_set = S.Reals)  # The solution is over real numbers
pprint(s1)
s2 = ConditionSet(x, x**2>4 , S.Reals)
pprint(s2)

⎧             2    ⎫
⎨x | x ∊ ℝ ∧ x  = 4⎬
⎩                  ⎭
⎧             2    ⎫
⎨x | x ∊ ℝ ∧ x  > 4⎬
⎩                  ⎭

19.4.7. Complement
The following script shows the use of Complement:-
This script is available on page 493 of the book

In [24]:
from sympy import *
s1 = FiniteSet(1,2,3,4)
s2 = FiniteSet(2,3)
s3 = Complement(s1, s2)
pprint(s3) # Output {1, 4}

{1, 4}

19.4.8. ImageSet (along with imageset function)
This topic is discussed in detail from page 493-495 of the book.
When a function say f is applied to a set say $A = {x| \ x ∈S}$ , then a new set say B is created such that $B = {f(x)| \ x ∈S}$. Then this set $B$ is the ImageSet of the set $A$. One advantage of using a set to represent the image is that it could represent those images which have infinite entries. For instance, equations like $sin(x) = 0$; $cos(x) = ½$ would have infinite solutions. Such images, that is, solutions can be easily represented by a set. The concept of ImageSet is explained very well in the SymPy docs as follows: “ Image of a set under a mathematical function. The transformation must be given as a Lambda function which has as many arguments as the elements of the set upon which it operates, for instance, 1 argument when acting on the set of integers or 2 arguments when acting on a complex region. This function is not normally called directly, but is called from imageset”.
This means that you can create the image of a set of values using either the ImageSet class or the imageset() function. The signature of ImageSet class is:

Init signature: ImageSet(Lamda, base_set)

Now the signature of imageset function is:

Signature: imageset(*args)
Docstring:Return an image of the set under transformation ``f``.
If this function can't compute the image, it returns an
unevaluated ImageSet object.
.. math::
    { f(x) | x \in self }

So the imageset() function takes a normal Python lambda function rather than a function returned by Lambda class.
This is shown in script below:-
This script is available on page 494 of the book

In [25]:
from sympy import imageset, Lambda, symbols, S
x,y = symbols('x y')
b = FiniteSet(1,2,3,4,8,16,36,49)# For intersection
# -------------Use sympy Lambda---------------
a = imageset(Lambda(x, 2*x), S.Integers)#Create imageset using Lambda from sympy
print('a->', a)
c = a.intersect(b)
print('c->', c)
print('c type->', type(c))
# ------------Use python lambda----------------
f = lambda x: x**2# Normal lambda of python
g = imageset(f, S.Integers) # Create imageset using normal python lambda
h = g.intersect(b)
print('g->', g)
print('h->', h)

a-> ImageSet(Lambda(x, 2*x), S.Integers)
c-> {2, 4, 8, 16, 36}
c type-> <class 'sympy.sets.sets.FiniteSet'>
g-> ImageSet(Lambda(x, x**2), S.Integers)
h-> {1, 4, 16, 36, 49}

19.4.9. Set operations in SymPy
Some common set functions in SymPy are shown in the following code:
This script is available on page 495 of the book

In [26]:
from sympy import imageset, Lambda, symbols, S
x,y = symbols('x y')
s1 = FiniteSet(1,2,3)
s2 = FiniteSet(3,4,5)
s3 = s1.union(s2)
print(s3)
s4 = s1 + s2 # Can use + operator for union
print(s4)
s5 = s1.intersection(s2)
print(s5)
s6 = s1-s2 # Can use - operator for difference
print(s6)
s7 = s1*s2 # Can use * operator for cartesian product
print(s7)
s8 = s1**2# Can use ** for cartesian product with itself
print(s8)
s9 = set(s1**2) # Expand the cartesian product to list all members
print(s9)

{1, 2, 3, 4, 5}
{1, 2, 3, 4, 5}
{3}
{1, 2}
{1, 2, 3} x {3, 4, 5}
{1, 2, 3} x {1, 2, 3}
{(1, 2), (3, 2), (1, 3), (3, 3), (3, 1), (2, 1), (2, 3), (2, 2), (1, 1)}

19.5. Matrices
SymPy provides extensive support for Matrices. It has a Matrix class which can take a list of lists. Each inner list is one row of the matrix. However, if you give only one list, that is, not a list of lists, then the list will be treated as a column matrix. This is clear from the following code:
This script is available on page 495 of the book

In [27]:
M1 = Matrix([[1,2,3], [4,5,6]]) # 2 x 3 matrix
display(M1)
M2 = Matrix([1,2,3,4]) #  4 x 1 Column matrix 
display(M2)

Matrix([
[1, 2, 3],
[4, 5, 6]])
Matrix([
[1],
[2],
[3],
[4]])

One can use Matrix.shape to get the shape of the matrix. You can get individual rows/ columns of the matrix using Matrix.row(row_num) or Matrix.column(col_num). You can also use row_del(row_num) or col_del(col_num) to delete rows/ columns. Similarly, you can use row_insert or col_insert for insertion of row/ column. You can also use operators like +/ - for addition/ subtraction. You can also use * for multiplication and you can raise a matric to power -1 to get inverse. You can also use Matrix.T to get transpose of a Matrix. All the above concepts are clarified in the following code:
This script is available on page 496 of the book

In [28]:
from sympy import *
from sympy.plotting import plot
from IPython.display import display
init_printing(use_latex='mathjax')
M1 = Matrix([[1,2,3], [4,5,6]]) # 2 x 3 matrix
print('M1 shape->', M1.shape) # Shape is 2 x 3
# M1 doesnt change on inserting row. New matrix created
M1_new = M1.row_insert(1, Matrix([[0,0,1]])) #Insert row [0,0,1]at index 1
print('M1_new ->', M1_new)
M2 = Matrix([[4,5,6], [7,8,9]]) # Another 2 x 3 matrix
M3 = M1 + M2
print('M1 + M2->', M3)
M4 = Matrix([[2,4], [3,5], [4,6]])# 3 x 2 matrix
M5 = M1*M4 # Multiply
print('M1*M4->', M5) # Note M5 is a square matrix so can find inverse

M6 = M5**(-1) # Get inverse of M5
print('Inverse of M5->', M6)
print('M5*M6->', M5*M6) #Confirms that M6 is inverse of M5
M7 = M1.T # M1 is 2 x 3. M7 is 3 x 2
print('Transpose of M1->', M7)

M1 shape-> (2, 3)
M1_new -> Matrix([[1, 2, 3], [0, 0, 1], [4, 5, 6]])
M1 + M2-> Matrix([[5, 7, 9], [11, 13, 15]])
M1*M4-> Matrix([[20, 32], [47, 77]])
Inverse of M5-> Matrix([[77/36, -8/9], [-47/36, 5/9]])
M5*M6-> Matrix([[1, 0], [0, 1]])
Transpose of M1-> Matrix([[1, 4], [2, 5], [3, 6]])

19.6. The Equality class and Eq
SymPy has a class Equality. This class can compare two objects. It has an alias Eq. So you may use Equality or Eq. The signature of Equality (or Eq) on Jupyter is as follows:

Init signature: Eq(lhs, rhs=0, **options)
Docstring:     
An equal relation between two objects. Represents that two objects are equal.  If they can be easily shown to be definitively equal (or unequal), this will reduce to True (or False).  Otherwise, the relation is maintained as an unevaluated Equality object.  Use the ``simplify`` function on this object for
more nontrivial evaluation of the equality relation. As usual, the keyword argument ``evaluate=False`` can be used to prevent any evaluation.
Difference between `solve()` and `solveset()`: This box shows the difference between use of `solve()` and `solveset(`) methods. (The method `solveset()` is discussed next, so you may read about the `solveset()` method before reading this box. Suppose you need to solve the equation $tan(x) = 0$. The method `solve()` gives only a single solution, that is, `[0]`, while `solveset(`) will give a solution of type: ${2.n.π | n ∊Z} U {2.n.π+ π| n∊Z$ The following code on Jupyter shows the difference:

This script is available in the box on page 497 of the book

In [29]:
from sympy import *
x = symbols('x')
# Equation to be solved is tan(x) = 0
#  Using solve()
s1 = solve(Eq(tan(x), 0), x)
pprint(s1)
# Using solveset()
s2 = solveset(Eq(tan(x), 0), x)
pprint(s2)

[0]
{2⋅n⋅π | n ∊ ℤ} ∪ {2⋅n⋅π + π | n ∊ ℤ}

19.7.1. solveset()
You can use solveset() to solve both equalities as well as inequalities. From the above signature of solveset(), it is clear that by default the domain over which the function is solved is complex. Following code shows use of solvset() to solve a quadratic equation of type $ax2 + bx + c = 0$ The following script shows this:-
This script is available on page 499 of the book

In [30]:
from sympy import *
a,b,c,x = symbols('a b c x')
q1 = a*(x**2) + 2*(b*x) + c
q2 = Eq(q1, rhs = 0)
s1 = solveset(q2, x)# Default for domain = S.Complexes
print('General solution->', s1)
q_sub = q1.subs([(a,1), (b, 2), (c, 5)])
s2 = solveset(q_sub, x, domain = S.Complexes)
print('Solution over Complex domain->', s2)
s3 = solveset(q_sub, x, domain = S.Reals)
print('Solution over Real domain->', s3)

General solution-> {-b/a - sqrt(-a*c + b**2)/a, -b/a + sqrt(-a*c + b**2)/a}
Solution over Complex domain-> {-2 - I, -2 + I}
Solution over Real domain-> EmptySet()

You can even plot a graph of function in SymPy because SymPy provides a plot() function. The details of the function are not discussed here (those interested may do ?plot on Jupyter). Basic signature of plot is as follows:

plot(expr, range, **kwargs)
``expr`` : Expression representing the function of single variable
``range``: (x, 0, 5), A 3-tuple denoting the range of the free variable.

The following script solves for roots of a cubic $x^3 – 4x^2 + x + 6 = 0$ and plots it:-
This script is available on page 499 of the book

In [31]:
% matplotlib inline
from sympy import *
from IPython.display import display
init_printing(use_latex='mathjax')
x = symbols('x')
exp1 = (x**3) - 4*(x**2) +x + 6
exp2 = Eq(exp1, rhs = 0)
sol1 = solveset(exp2, x)
plot(exp1, (x,-2, 4 ))
Out[31]:
<sympy.plotting.plot.Plot at 0x5f14c90>

In fact solveset() function even express those solutions which may be infinite. For instance, see the code below for $sin(x) =0$.
This script is available on page 499 of the book

In [32]:
exp1 = sin(x)
exp2 = 0
s1 = solveset(Eq(exp1, exp2), x)
display(s1)
$$\left\{2 n \pi\; |\; n \in \mathbb{Z}\right\} \cup \left\{2 n \pi + \pi\; |\; n \in \mathbb{Z}\right\}$$

Note that solveset() returns a set object and so it takes care of different types of output. For instance, if there is no solution, then the set returned will be an empty set. The solveset() method provides the solution in the form of a set. There is also a solvify() method to provide solutions in form of a list. This is shown in script below:-
This script is available on page 500 of the book

In [33]:
from sympy import *
from sympy.solvers.solveset import solvify# Have to import solvify
from IPython.display import display
init_printing(use_latex='mathjax')
x = symbols('x')
exp1 = (x**3) - 4*(x**2) +x + 6
exp2 = Eq(exp1, rhs = 0)
sol1 = solvify(exp2, x, domain = S.Complexes)
print(sol1)

[-1, 2, 3]
Advanced: (For details see the corresponding box on page 500 in the book) The distinction between a simple `FiniteSet` and a `FiniteSet of ordered tuples` is subtle but very important. This is because in general a `FiniteSet` is “`unordered`” so a `FiniteSet` say `{1, 2, 3, 4}` is same as `{4, 3, 2, 1}`. However, a `FiniteSet` of ordered tuple is of type `{(1, 2, 3, 4)}` and is “`ordered`” because inside the `FiniteSet`, there is a tuple. So you can “map” the “order of variables” to the “order in the solution”. For instance, you may decide in a three-dimensional problem that the order is `{(x, y, z)}`. Therefore, if you get say `{(1, 2, 3)}`, then you have `x =1, y =2 and z = 3`. This is clear from the following code:

This script is available in the box on page 501 of the book

In [34]:
from sympy import *
x,y,z = symbols('x y z')
s1 = FiniteSet(x,y,z) # Unordered
s2 = FiniteSet(z,x,y) # Unordered
s3 = FiniteSet((x,y,z)) # Ordered
s4 = FiniteSet((z,x,y)) # Ordered
s5 = s1 -s2 # Gives EmptySet()
print(s5)
s6 = s3 - s4 # Gives {(x, y, z)} \ {(z, x, y)}
print(s6)

EmptySet()
{(x, y, z)} \ {(z, x, y)}

19.8. The linsolve() method
(For detailed explanation, see the book. Here only the script is given)

Consider the following system of three linear equations in three variables: x, y and z $x + y + z = 5$,
$2x + 3y + 5z = 8$ and
$4x + 5z = 2$.
First method:
In the first method the parameters are provided to the linsolve() method in the form of (A,b) where A and b are matrix. This is shown in the following code:
This script is available on page 502 of the book

In [35]:
from sympy import Matrix, S, linsolve, symbols
#x + y + z = 5, 2x + 3y + 5z = 8, 4x + 5z = 2
x, y, z = symbols("x, y, z")
A = Matrix([[1, 1, 1], [2, 3, 5], [4, 0, 5]])
b = Matrix([5, 8, 2])
sol = linsolve((A, b), [x, y, z])
print(sol)# Solution is {(3, 4, -2)}

{(3, 4, -2)}

Second method:
Here the three linear equations are provided as comma separated values to the system parameter of the linsolve() method as shown in the following code:
This script is available on page 503 of the book

In [36]:
from sympy import *
from sympy.solvers.solveset import linsolve
x, y, z = symbols('x y z')
my_sol = linsolve([x + y + z - 5, 
                   2*x + 3*y + 5*z - 8, 
                   4*x + 5*z - 2], (x, y, z))
print(my_sol)

{(3, 4, -2)}

19.9.1. Differentiation
For differentiation, SymPy has a function diff(expression, variable_1, variable_2…). Here expression is the expression to be differentiated and variable_1, variable_2 are the variables with respect to which the differentiation is to take place. So you can do any number of differentiations at once. This will become clear from the following example:
This script is available on page 503 of the book

In [37]:
from sympy import *
x, y, z = symbols('x y z')
exp1 = x*((x**3) + (y**4) + (z**5))
exp2 = diff(exp1, x)
print(exp2) # exp2-> 4*x**3 + y**4 + z**5
exp3 = diff(exp1, x, y)
print(exp3) # exp3-> 4*y**3
exp4 = diff(exp1, x, y, y)
print(exp4) # exp4-> 12*y**2

4*x**3 + y**4 + z**5
4*y**3
12*y**2

You can differentiate an expression with respect to another expression also. See the following code:

In [38]:
from sympy import Symbol
x = Symbol('x')
print((sin(x)**2).diff(x))# Differentiate wrt x
print((sin(x)**2).diff(sin(x)))# Differentiate wrt sin(x)

2*sin(x)*cos(x)
2*sin(x)

SymPy also provides a Derivative class for derivation.
This script is available on page 504 of the book

In [39]:
from sympy import *
x, y = symbols('x y')
exp1 = 2*(x**2) + 5*x*y
sol = Derivative(exp1, x)
print(sol)

Derivative(2*x**2 + 5*x*y, x)

19.9.2. Integration
Similarly to integrate you can use the function integrate(expression, var_1, var_2…). This will be clear from the following code:
This script is available on page 504 of the book

In [40]:
from sympy import *
x, y, z = symbols('x y z')
exp1 = sin(x) + log(y) + z
exp2 = integrate(exp1, x)
print(exp2) # exp2-> x*z + x*log(y) - cos(x)
exp3 = integrate(exp1, x, y)
print(exp3) # exp3-> x*y*log(y) + y*(x*z - x - cos(x))

x*z + x*log(y) - cos(x)
x*y*log(y) + y*(x*z - x - cos(x))

You can also use the integrate function for integrating definite integrals. The format is integrate(expression, (var_1, LowLmt, UppLmt), (var_2, LowLmt, UppLmt)…..).
Note: ∞ (infinity) in SymPy is oo (the lowercase letter “o” twice). You can use SymPy for plotting a function also. This will be clear from following code:
This script is available on page 504 of the book

In [41]:
from sympy import symbols
from sympy.plotting import plot
from IPython.display import display
init_printing(use_latex='mathjax')
%matplotlib inline
exp1 = integrate(x**x, (x, 0, 1))
display(exp1)
plot(x**x, (x,0,1))
$$\int_{0}^{1} x^{x}\, dx$$
Out[41]:
<sympy.plotting.plot.Plot at 0x302c0f0>

19.9.3 Finding limits (Lim)
SymPy can also be used to get limits of functions. You can solve an equation of type $ \lim_{x \to 0} f(x)$
For instance, it is well known that
$ \lim_{x \to 0} (x + \frac{1}{x}) = e$
You can solve this on SymPy as follows:
This script is available on page 505 of the book

In [42]:
x = symbols('x')
exp1 = (1 + 1/x)**x # Function whose limit to be found
display('f(x)->', exp1)
exp2 = limit(exp1, x,oo ) # limit at infinity
display('lim f(x) ->',exp2)

'f(x)->'
$$\left(1 + \frac{1}{x}\right)^{x}$$ 
'lim f(x) ->'
$$e$$

In SymPy you can evaluate limit from both sides, that is, positive side as well as negative side. Following example clarifies it:
This script is available on page 505 of the book

In [43]:
exp1 = (1/(x-3)) # Function whose limit to be found
exp2 = limit(exp1, x,3, '+' ) # limit from positive direction
display('lim f(x) ->',exp2) # Gives oo (+ infinity)
exp3 = limit(exp1, x,3, '-' ) # limit from negative direction
display('lim f(x) ->',exp3) # Gives -oo (- infinity)

'lim f(x) ->'
$$\infty$$ 
'lim f(x) ->'
$$-\infty$$

19.9.4. Ordinary Differential equations (ODE)
SymPy provides a dsolve() method for solving differential equations. This method is part of the ODE module. The dsolve() method has many parameters as can be seen by its signature. Here it is not possible (nor necessary) to discuss all the details. The complete signature of dsolve() is:

dsolve(eq, func=None, hint='default', simplify=True, ics=None, xi=None, eta=None, x0=0, n=6, **kwargs)

19.9.5. Solving ODE for an undamped and damped harmonic oscillator
(For detailed explanation see the book)
$m*ÿ(t) + k*y(t) = 0$ # Undamped oscillator ie c = 0
$m*ÿ(t) +c* ẏ(t) + k*y(t) = 0$ # Damped oscillator
Where

  • $ÿ(t) = d^2y(t)/ dt^2$
  • $ẏ(t) = dy(t)/dt$
  • $t = time$
  • $m = mass$
  • $k = spring constant$
  • $c = friction constant$ For those unfamiliar with oscillators, the terminology doesn’t matter. Just think of the above two differential equations. The following script provides the solution to these ODE:
    This script is available on page 507 of the book
In [44]:
from sympy.solvers.ode import dsolve
from sympy import *
t, m, k, c = symbols('t m k c')
y = Function('y')(t)
y_ = Derivative(y, t) # y_ is first derivative of y(t) wrt to t
y__ = Derivative(y_, t) # y__ is second derivative of y(t) wrt to t
# Undamped oscillator m is mass k is spring constant
eq_ud = m*y__ + k*y
sol_ud = dsolve(eq_ud)
print('Undamped oscillator->', sol_ud)
print('type of differential eq->', classify_ode(eq_ud))  # The classify_ode() method gives "type" of ODE
# Damped Oscillator force of friction Ff = -cy_
eq_dd = m*y__ + c*y_ + k*y
sol_dd = dsolve(eq_dd)
print('Damped oscillator->', sol_dd)
print('type of differential eq->', classify_ode(eq_dd))

Undamped oscillator-> Eq(y(t), C1*exp(-t*sqrt(-k/m)) + C2*exp(t*sqrt(-k/m)))
type of differential eq-> ('nth_linear_constant_coeff_homogeneous', '2nd_power_series_ordinary')
Damped oscillator-> Eq(y(t), C1*exp(t*(-c - sqrt(c**2 - 4*k*m))/(2*m)) + C2*exp(t*(-c + sqrt(c**2 - 4*k*m))/(2*m)))
type of differential eq-> ('nth_linear_constant_coeff_homogeneous', '2nd_power_series_ordinary')

19.11 Exercise
The equation for radioactive decay is given by:
$dn(t)/dt= -λ\times t$ Set up the above equation in SymPy and also find its solution.
Hint:
This script is available on page 509 of the book

In [45]:
from sympy import *
init_printing()
# You can use 'lambda' to print the Greek letter for lambda
# Use capital L since small l is difficult to read
L = symbols('lambda')
# t is for time
t = symbols('t')
# n(t) is the number of particles at time t.
n = Function('n')(t)
# n_ means differential of n with respect to t ie dn(t)/dt
n_ = n.diff(t)
eq_decay = Eq(n_, -L*n)
# Print the equation for radioactive decay
pprint(eq_decay)
# Use dsolve to get solution of the ode
decay_sol = dsolve(eq_decay)
pprint(decay_sol)

d                 
──(n(t)) = -λ⋅n(t)
dt                
           -λ⋅t
n(t) = C₁⋅ℯ

Exercise (b)

Using SymPy, solve the following ODE:

$\frac{dy}{dx} = \frac{(1 - y)^{4/3}}{y}$
Hint:
This script is available on page 509 of the book

In [46]:
from sympy import symbols, Eq, Derivative, Rational, dsolve, pprint
y = symbols('y')
x = symbols('x')
eq = Eq(Derivative(y(x),x), (1 - y(x))**Rational(4,3) / y(x))
pprint(dsolve(eq))

y(x)                               
 ⌠                                 
 ⎮           y                     
 ⎮   ────────────────── dy = C₁ - x
 ⎮   3 ________                    
 ⎮   ╲╱ -y + 1 ⋅(y - 1)            
 ⌡

Beyond text book
This topic is given on page 509 of the book
(b) Getting to know the type of ODE and solving Bernoulli ODE
The ODE module of SymPy has a method classify_ode(). The signature of the method is:

sympy.solvers.ode.classify_ode(eq, func=None, dict=False, ics=None, **kwargs)
Returns a tuple of possible dsolve() classifications for an ODE.

A Bernoulli ODE is of type:
$ \frac{dy}{dx} + P(x).y = Q(x).y^n $
The following code uses SymPy to solve the following two Bernoulli ODEs:
$ eq1 = \frac{dy}{dx} + xy - xy^3 $
$ eq2 = 6\frac{dy}{dx} -2y + xy^4$
The following script classifies the two ODEs and also solves them:
This script is available on page 511 of the book

In [47]:
from sympy.solvers.ode import dsolve, classify_ode
from sympy import *
init_printing(use_latex='mathjax')

x = symbols('x')
y = Function('y')(x)
y_ = Derivative(y, x)
# Both eq1 and eq2 are Bernoulli ODE
eq1 = y_ + x*y + x*y**3
eq2 = 6*y_ -2*y -x*y**4
# The classify_ode() method gives a tuple of possible classification of ODE
eq1_class = classify_ode(eq1)
eq2_class = classify_ode(eq2)
print(eq1_class)
print(eq2_class)
sol_ode1 = dsolve(eq1)
sol_ode2 = dsolve(eq2)
pprint(sol_ode1)
pprint(sol_ode2)

('separable', '1st_exact', 'Bernoulli', '1st_power_series', 'lie_group', 'separable_Integral', '1st_exact_Integral', 'Bernoulli_Integral')
('Bernoulli', '1st_power_series', 'lie_group', 'Bernoulli_Integral')
⎡              ____________               ____________⎤
⎢             ╱    -C₁                   ╱    -C₁     ⎥
⎢y(x) = -    ╱  ────────── , y(x) =     ╱  ────────── ⎥
⎢           ╱         ⎛ 2⎞             ╱         ⎛ 2⎞ ⎥
⎢          ╱          ⎝x ⎠            ╱          ⎝x ⎠ ⎥
⎣        ╲╱     C₁ - ℯ              ╲╱     C₁ - ℯ     ⎦
                     1              
y(x) = ─────────────────────────────
            ________________________
           ╱ ⎛               x⎞     
          ╱  ⎜     (-x + 1)⋅ℯ ⎟  -x 
       3 ╱   ⎜C₁ + ───────────⎟⋅ℯ   
       ╲╱    ⎝          2     ⎠

Beyond text
c. Using sample ODE (with solutions) from file test_ode.py
This topic is given on page 511 of the book
SymPy has a file test_ode.py which provides many ODE and their solutions and compares the ODE to its solution. This is a huge file with about 3000 lines of code. One need not go into the details of the file. This file however provides a large number of ODE and their solutions and compares them. To do this the file has a large number of functions. The following script considers the simplest one, that is:

def test_linear_2eq_order1():
    x, y, z = symbols('x, y, z', cls=Function)
    k, l, m, n = symbols('k, l, m, n', Integer=True)
    t = Symbol('t')
    x0, y0 = symbols('x0, y0', cls=Function)
    eq1 = (Eq(diff(x(t),t), 9*y(t)), Eq(diff(y(t),t), 12*x(t)))
    sol1 = [Eq(x(t), 9*C1*exp(6*sqrt(3)*t) + 9*C2*exp(-6*sqrt(3)*t)), \
    Eq(y(t), 6*sqrt(3)*C1*exp(6*sqrt(3)*t) - 6*sqrt(3)*C2*exp(-6*sqrt(3)*t))]
    assert checksysodesol(eq1, sol1) == (True, [0, 0])

In the above code, the ODE is given in line number 34 and its solution is given in lines 35–36. The name of the function, that is, test_linear_2eq_order1() shows that this particular function is for a system of two linear differential equations of order 1 (The rest of the file contains many other cases).
The following script uses the above differential equation and sees if it actually gives the indicated result:
This script is available on page 512 of the book

In [48]:
from sympy.solvers.ode import dsolve, checksysodesol
from sympy import *
#1---------------define symbols and functions----------------
# Following symbols and Function copied from test_ode.py
C0, C1, C2, C3, C4, C5, C6, C7, C8, C9, C10 = symbols('C0:11')
x, y, z = symbols('x:z', real=True)
f = Function('f')
g = Function('g')
h = Function('h')
x, y, z = symbols('x, y, z', cls=Function)
k, l, m, n = symbols('k, l, m, n', Integer=True)
t = Symbol('t')
# 2---------------copy the ode from test_ode.py
eq1 = (Eq(diff(x(t),t), 9*y(t)), Eq(diff(y(t),t), 12*x(t)))
# 3------- Solve the ode using dsolve()------------
sol2 = dsolve(eq1)
print(sol2)

[Eq(x(t), 9*C1*exp(-6*sqrt(3)*t) + 9*C2*exp(6*sqrt(3)*t)), Eq(y(t), -6*sqrt(3)*C1*exp(-6*sqrt(3)*t) + 6*sqrt(3)*C2*exp(6*sqrt(3)*t))]